function output = rhythm_gpu(y);
if nargin < 1, y = wavread('technobeat.wav'); end

bandlimits=gsingle([0 200 400 800 1600 3200]);
numBands = length(bandlimits);
maxSigFreq=44100;
y = gsingle(y);
hannWinLength = .4;
hannlen = hannWinLength*2*maxSigFreq;
yFFTed = fft(y);
lenY = length(yFFTed);
bandRanges = gzeros(12,1);

% Using band limits to find the ranges
i = [];
gfor i = 1:numBands-1
    bandRanges(2*(i)-1) = floor(bandlimits(i)/maxSigFreq*lenY/2)+1;
    bandRanges(2*i) = floor(bandlimits(i+1)/maxSigFreq*lenY/2);
    
gend
% The end cases of the band limits are special cases
bandRanges(2*numBands-1) = floor(bandlimits(numBands)/maxSigFreq*lenY/2)+1;
bandRanges(numBands*2) = floor(lenY/2);

output = gzeros(lenY,numBands); %for storing the output
% Using the frequency bands to format the output
for i = 1:numBands
    start = bandRanges(2*i-1);
    stop = bandRanges(2*i);
    output(start:stop,i) = yFFTed(start:stop);
    output(lenY+1-stop:lenY+1-start,i) = yFFTed(lenY+1-stop:lenY+1-start);
end
output(1,1)=0;

% signal in time domain
sigTime = gones(size(output)); 
% signal in frequency domain
sigFreq = gones(size(sigTime)); 
% for the hann 
hann = zeros(lenY,1); 

%for the (relatively) small value of hannlen, it's slower to use the GPU
for a = 1:hannlen
    hann(a) = (cos(a*pi/hannlen/2)).^2;
end
%now converting to GPU
hann = gsingle(hann);

% Going from frequency domain and then back with absolute values in the time domain;
gfor k = 1:numBands
    sigFreq(:,k) = fft(abs(real(ifft(output(:,k)))));
gend

%sigFreq = fft(abs(sigTime));

% Half-Hanning FFT * Signal FFT. And then back to time domain
i = [];
gfor i = 1:numBands
    output(:,i) = real(ifft(sigFreq(:,i).*fft(hann)));
gend
